####复杂热图###########################
rm(list = ls())
options(stringsAsFactors = F)

library(clusterProfiler)
library(msigdbr)
library(GSVA) 
library(GSEABase)
library(pheatmap)
library(limma)
library(colorRamps)
library(ComplexHeatmap)
library(vegan)
library(circlize)


setwd('父级路径')

##提取自己的样本 下面gcr 开头是我写的理解 
FPKM=read.table("convert_exp.txt",sep="\t",header=T,check.names=F, row.names = 1)
####  GSVA  ####
dat <- as.matrix(log2(FPKM+1))
##这一步是提取自己的基因集
genelist <- read.table("工作簿3.csv",sep=",",header=T,check.names=F)
geneset <- lapply(genelist, function(col) col[!is.na(col) & col != ""])

##linux系统
gsva_mat <- gsva(expr=dat, 
                 gset.idx.list=geneset, 
                 kcdf="Gaussian" ,#"Gaussian" for logCPM,logRPKM,logTPM, "Poisson" for counts
                 verbose=T, 
                 parallel.sz = parallel::detectCores())#调用所有核
##win系统
counts=read.table("GSE126044_counts.txt",sep="\t",header=T,check.names=F, row.names = 1)
####  GSVA  ####
dat <- as.matrix(counts)
genelist <- read.table("工作簿3.csv",sep=",",header=T,check.names=F)
geneset <- lapply(genelist, function(col) col[!is.na(col) & col != ""])

##linux系统
gsva_mat2 <- gsva(expr=dat,
                  gset.idx.list=geneset,
                  kcdf="Gaussian" ,#"Gaussian" for logCPM,logRPKM,logTPM, "Poisson" for counts
                  verbose=T,
                  parallel.sz = parallel::detectCores())#调用所有核
##win系统
merge <- cbind(gsva_mat2,gsva_mat)
#gcr 刚才说这里是错误的需要用合并后的数据 
write.csv(merge,"gsva_go_matrix.csv")


# 生成图
# gcr  ConsensusClusterResults.k=3.consensusClass.csv 这个文件是用于聚类使用的吧
cluser = read.csv(file = 'ConsensusClusterResults.k=3.consensusClass.csv', sep = ',', header = F)
rt <- read.csv(file = 'gsva_go_matrix.csv', sep = ',', header = T, row.names = 1)
rownames(cluser) <- cluser[,1]
names(cluser)[2] <- "Cluster"
cluster1_rows <- rownames(cluser[cluser$Cluster == "1", ])#看一下1有多少，这里是221
rt_intersect1 <- rt[, intersect(colnames(rt), cluster1_rows)]
##再来一遍
cluster2_rows <- rownames(cluser[cluser$Cluster == "2", ])#看一下2有多少，这里是276
rt_intersect2 <- rt[, intersect(colnames(rt), cluster2_rows)]
cluster3_rows <- rownames(cluser[cluser$Cluster == "3", ])#看一下3有多少，这里是176
rt_intersect3 <- rt[, intersect(colnames(rt), cluster3_rows)]

merged <- cbind(rt_intersect1, rt_intersect2, rt_intersect3)
merged<-as.matrix(merged)
merged <- decostand(merged,"standardize",MARGIN = 1)#标准化
# gcr  后期这里面颜色采用前端调色板输入
col_fun = colorRamp2(c(-2, 0, 2), c("#1f60a4","white","#a40d20"))#,space = "RGB"
col_Cluster1 <- c("1" = "#B49FDA", "2" = "#72AFD3", "3" = "#5FDFFF")
cluser <- cluser[colnames(merged),]
column_ann <- HeatmapAnnotation(Cluster=cluser$Cluster,
                                col = list(Cluster = col_Cluster1)
)

# draw(column_ann)
# gcr 上面这个我注释掉了 下面的这个方法我采用前端动态参数输入 
# gcr 目前也就是两个输入文件对吧 就是上面的 convert_exp.txt (自己样本 ) 工作簿3.csv 这两个文件是输入文件 生成的gsva_go_matrix.csv 文件当做 下面画图的输入 对吧
p1 <- Heatmap(merged,
              col = col_fun,
              cluster_rows = T,##打开行聚类
              cluster_columns = T,##关闭列聚类
              show_row_dend = T, ## 显示行树状图
              show_column_dend = T,## 显示列树状图
              show_row_names = T,##行名
              show_column_names = T, ##列名
              #column_names_rot = 2,##列名旋转方向
              row_names_rot = 0,##行名旋转方向角度 0 默认
              column_names_side = "top",#列名位置在顶部还是底部
              column_dend_side = "top",#柱状图应该放在热图的顶部还是底部
           
              top_annotation = column_ann
)

p1
# 复杂热图
Cairo::CairoTIFF(file="p1.tiff", width=8, height=8,units="in",dpi=150)
p1
dev.off()

pdf("Heatmap.pdf",width = 8,height = 8)
p1
dev.off()

